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NUMERICAL COMPUTATION OF TWO-DIMENSIONAL VISCOUS 
BLUNT BODY FLOWS WITH AN IMPINGING SHOCK 


John C. Tannehill and Terry L. Holst 
Iowa State University 


SUMMARY 


Two-dimensional viscous blunt body flows with an impinging shock 
have been computed using a "time-dependent" finite-difference method which 
solves the complete set of Navier-Stokes equations for a compressible 
flow. For low Reynolds number flows, the entire flow field, including 
the bow shock and impinging shock, has been "captured" in the computa- 
tion. For higher Reynolds number flows, the bow shock is treated as a 
discontinuity across which the Rankine-Hugoniot equations are applied, 
while the boundary layer and interaction regions are "captured" as 
before. Using this latter "shock-fitting" approach, a Type III shock 
interaction flow field has been computed with flow conditions 
corresponding to the Space Shuttle Orbiter freestream conditions at 
61 km (200,000 ft). 
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NOTATION 

specific heat at constant pressure 
constant in Sutherland's equation 
diameter of cylinder 
specific internal energy 
total energy 

unit vector in radial direction 
unit vector in transverse direction 
coefficient of thermal conductivity 
local radius of curvature 
Knudsen number 
Mach number 

outward unit normal to bow shock 
pressure 

Prandtl number, c u/k 
P 

heat flux vector 

local. radius of bow shock 

local radial shock velocity 

dr /ae 

s 

gas constant 

Reynolds number based on diameter 
time 

transformed time 

tangential velocity component 
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U = local shock velocity 

s 

v = normal velocity component 

= fluid velocity normal to shock 
y = transformed coordinate given by Eq. (16) 

y = transformed coordinate given by Eq. (14) 

z = transformed coordinate given by Eq. (16) 

z = transformed coordinate given by Eq. (14) 

Oi = transformation parameter in Eq. (Al) 

P = stretching factor in Eq. (16) 

j3 = stretching factor in Eq. (Al) 

■y = specific heat ratio 

6 = local distance between body and outer boundary 

T) = coordinate measured normal to body 

8 = angle measured from stagnation streamline 

fi = viscosity coefficient 

5 = coordinate measured along body 

p = density 

T. . ..= shear stress tensor 

ij 




,,aU' e 
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INTRODUCTION 


An extraneous shock impinging on a blunt body in a hypersonic flow 
has been observed to greatly increase both the heat transfer rate and 
pressure near the impingement point. In fact, Hains and Keyes^ have 
measured peak heating rates up to 17 times the ordinary stagnation 
point rate and pressure peaks up to 8 times the freestream pitot pressure 
as a result of shock impingement. Flow fields of this type will occur 
on the Space Shuttle and other maneuverable re-entry vehicles. For 
example, high heating rates and pressures can be expected on the Space 
Shuttle Orbiter at the point where the bow shock from the nose intersects 
the blunt leading edge of the wing if the sweep angle is not large 
enough. 

The intense heating and high pressures described above occur over 
a small region where a disturbance, originating at the intersection of 
the impinging shock and bow shock, strikes the body. The disturbance 

may be a free shear layer, a supersonic jet, or a shock depending on the 

_ , 2 
location of the impinging shock and the shape of the body. Edney has 

described six different types of shock interference patterns, which 

include the various kinds of disturbances. The shock interference 

pattern which produces the maximum heating rates and pressures is Type IV 

which is shown in Fig. 1. In this type of interference pattern, the 

disturbance is a supersonic jet which is embedded in the subsonic 

portion of the flow field. 

During the present study,' two-dimensional shock impingement flow 
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fields have been computed using a "time-dependent," finite-difference 
method which solves the complete set of Navier-Stokes equations for a 
compressible flow. The major advantage of the "time-dependent" method is 



Fig. 1. Type TV shock interference pattern. 

that the resulting unsteady Navier-Stokes equations are a mixed set of 
hyperbolic-parabolic equations for both subsonic and supersonic flows. 

As a result, a very complicated flow field, such as the one shown in 
Fig. 1 where both subsonic and supersonic regions are present, can be 
calculated as an initial -value problem. An additional advantage is that 
since the Navier-Stokes equations are solved in a conservative manner, 
shocks are automatically allowed to form without previous knowledge of 


their location or even existence. 
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For low Reynolds number flows, the entire flow field, including 
the bow shock, can be "captured" using the so-called "shock-capturing" 
approach. The computational domain for this type of calculation is shown 
in Fig. 2. For higher Reynolds number flows, it is not practical to 



Fig. 2. Computational domain for "Shock-Capturing" method. 

"capture" the bow shock because of the numerical difficulties associated 
with the large gradients at the bow shock. Instead, it is more 
convenient to treat the bow shock as a discontinuity, across which the 
Rankine-Hugoniot equations can be applied, while leaving the boundary 
layer and interaction regions to be captured as before. This latter 
approach is the so-called "shock-fitting" method. The computational 
domain for this method is shown in Fig. 3. In the present study, both 
the "shock-capturing" and "shock-fitting" methods have been used to 
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COMPUTATIONAL 



Fig. 3. Computational domain for "Shock- Fit ting" method . 

compute two-dimensional viscous blunt body flows with an impinging shock. 
Thus far, numerical difficulties have prevented the computation of 
Edney's Type IV shock interference patterns. . However, Type III inter- 
actions have been successfully computed in this study. 
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GOVERNING EQUATIONS 

The fundamental governing equations for an unsteady flow without 

body forces or external heat additions can be written in conservation 

law form for a two-dimensional body, intrinsic coordinate system (see 
3 

Fig. 4) as 


st + sl + d?i +H_0 


( 1 ) 


where 


U = (1 + KT|) 


P 

pu 

pv 

E 


( 2 ) 


F = 


pu 


P + PU - T^r 
puv - T_ 

Eu + pu + q- - u T 


§1 " vT sri 


( 3 ) 


G = (1 + KID 


pv 

puv - T 


r\ 5 


O + pv - T, 


rm 


Ev + pv + q. - uT - - vt 

H n td s nn 


( 4 ) 



6 


H = 


0 

K(puv - V 

-K(p + pu 2 - T ?g ) 
0 


(5) 


and E is the total energy given by 


E = 



+ 



( 6 ) 


In addition to the above conservation equations, an equation of state in 
the form 


P = p(e,P) (7) 

must be specified. For a perfect gas, this equation can be written as 

P = (7 - 1) Pe (8) 

For the case of air in chemical equilibrium, approximate curve fits 
are available for Eq. (7) in Refs. 4 and 5. 



Fig. 4. Two-dimensional body intrinsic coordinate system. 
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In the present study, the Navier-Stokes expressions were used for 
the shearing stress tensor and heat flux vector; thus the components of 
the shearing stress tensor and heat flux vector are 


T SS “ - 3 11(8 1? + e W 

t tfi' “ e rrn ■ 3 H<e SS + 


t 5>i |je ?T1 

where 



( 1 \ /Sv\, /_Ku \ 

e ?ri ^ 1 + kti ) l S§ I Sri l 1 + KTiy 

and 

-k Sr 

_ 1 + KTI 

. Sr 

q n " " k an 


(10) 


(ii) 


In order to complete the system of equations, it is necessary to 
specify expressions for the viscosity (p,) and the coefficient of thermal 
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conductivity (k) . For all of the present computations, Sutherland's 
equation 


H = 



1 + C/T 
1 + C/T 


( 12 ) 


was used to determine the viscosity, and the coefficient of thermal 
conductivity was computed by assuming a constant Prandtl number: 


k = 


ULi± 

(y - l)Pr 


(13) 


Approximate curve fits for p, and k suitable for "time-dependent" computa- 
tions are presently being developed for the case of air in chemical 
equilibrium. 

Two independent variable transformations are applied to the governing 
equations listed above. The first transformation maps the computational 
domain shown in Fig. 4 into a rectangular region in the transformed 
(y, z) plane. The outer computational boundary in Fig. 4 is either the 
bow shock (for a "shock-fitting" computation) or a freestream boundary 
(for a "shock-capturing" computation). The equations of the first 
independent variable transformation are: 


y = 5 

Z = 1 - J- (14) 

t = t 


where 6 = 6(5, t) for a "shock- fitting" computation and 6 = 6(5) for a 
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"shock-capturing" computation. The relations between the partial derivatives 
are given by 


d_ 


= -^ + (1 - z) 
dy 



d_ 

dri 


1 

6 


d_ 

dz 


d_ 

at 




( 15 ) 


The second independent variable transformation stretches the computa- 
tional mesh in the direction normal, to the body so that the mesh may be 
refined either near the body or near the outer computational boundary. 

By refining the mesh near the body, it is possible to describe the 
boundary layer more accurately. On the other hand, it may be desirable 
to refine the mesh near the outer boundary for shock impingement computa- 
tions, although the results from initial computations tend to dispute 
this. The equations of this transformation, which were previously used 
by Li^, are 

y = y 



t = t 


The amount and type of stretching is determined by the value of P. If P 
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equals zero there is no stretching. If 3 is greater than zero the mesh 
is refined near the body, and if (3 is less than zero the mesh is refined 
near the outer computational boundary. 

The relations between the partial derivatives for this second trans- 
formation are given by 



d_ 

dz 



(17) 




An alternative to the present transformation is given in Appendix A. 

This alternate transformation makes use of a logarithmic function which 
has been found to give better computational results when an impinging 
shock is present. 

After employing the two transformations given by Eqs. (14) and (16), 
the final computational grid in the (y,z) plane is shown in Fig. 5 and 
the corresponding grid in the physical plane (§, r l) is shown in Fig. 6 
for P > 0. The final forms of the governing equations are 


dU 

at 


dF 


dy dz 


+ ^+ ^-(-H = 0 


(18) 
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Fig. 6. Physical plane ( p > 0) 
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U = [6 + 6 2 k(L - z)] 
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NUMERICAL SOLUTION OF EQUATIONS 


Finite-Difference Scheme 

7 8 

MacCormack's finite-difference scheme * is used to solve the govern 
ing equations at each interior grid point. This explicit, two-step 
scheme has second-order accuracy in both space and time. When 
MacCormack’s algorithm is applied to Eq. (18), the following predictor- 
corrector equations result 


u .*: 1 . u " . 

J,k J,k 


^ (F n 
Ay ^j.-k+l 


F n 


.) - 


— fG n * - 
Az CG j+l,k 


G ",k> - 4t H j,k < 27 > 


and 


u n+ i - i u" .4. o.t 1 - 
j,k 2 j,k j,k 


?. n + l A - G. n f. N ) - At 

\ j. k J,k-1 J Az\^j,k J-l,k/ 




(28) 

where z = jAz, y = kAy, t = nAt and F*! , = FCU 1 ? , ), F™ + ? = F(U. n ^), etc. 

j,K j,k j,k j,k 

Note that the spatial derivatives in the predictor step are approximated 

by forward difference, while in the corrector step they are approximated 

by backward differences. The shear stress and heat flux terms appearing 

in F, G, and H are evaluated using backward differences in the predictor 

step and forward differences in the corrector step. The net result is 

a central difference approximation for the shear stress and heat flux 
9 

terms . 


Using the above finite-difference equations, the computation is 
advanced in time from the initial conditions until the "steady-state' 
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solution is reached. After each computational step, the flow variables 
are obtained at each interior grid point from the U vector 


- 

U 1 


P 

U 2 

= [6 + 6 2 K(1-z)] 

pu 

U 3 


pv 



E 

U 4 

— 


L_ 


(29) 


in the following manner: 


P - flj/ [6 + 6 2 K(1 -z)] 


u - U 2 /U t 


v - 


e - U./u. 

4 1 


2 2 
u + v 


P = p(e,p) 


(30) 


T = T(e, p) 

During the present study, the latter two expressions were evaluated 
using perfect gas relations. However, for the case of air in chemical 
equilibrium, they could easily be evaluated using the approximate curve 
fits appearing in Refs. 4 and 5. 


Boundary Conditions 


The flow conditions at the outer computational boundary in Fig. 7 



16 


OUTER 



BOUNDARY 

Fig. 7 Computational boundaries 

are necessarily different for the "shock-capturing" and "shock-fitting" 
methods. For the "shock-capturing" method, freestream conditions are 
maintained at all grid points along the outer computational boundary 
below the impingement point. For grid points above the impingement point, 
the flow variables are set equal to the conditions which exist behind an 
oblique shock at the desired shock impingement angle. The impingement 
point is placed, for simplicity, halfway between any two adjacent grid 
points along the outer computational boundary. 

When the "shock- fitting" method is employed, the flow conditions 
at the outer computational boundary are those conditions which exist 
immediately downstream of the bow shock as determined by the Rankine- 
Hugoniot relations. Consequently, it is necessary to include logic 
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which will permit this boundary to move with the bow shock as the latter 
moves toward its "steady-state" position. The approach used here is 
somewhat similar to the approach previously used by Thomas, et al.^ 
and Kutler, et al.^ in their inviscid steady flow computations. Their 
predictor-corrector approach has been modified for the present unsteady 
computations . 

The coordinate system used for the "shock-fitting" procedure 
is shown in Fig. 8. The local velocity of the shock is given 



Fig. 8 Notation for "shock- fitting" procedure. 


by 

U = U n (31) 

s s s 

where n^ denotes the outward unit normal to the shock given by 
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i - ( r /r j i Q 

r \ s e s / 6 

[1 + (r /r ) 2 ] % 

s 8 s 


(32) 


The magnitude of the local shock velocity can be related to the radial 
shock velocity (r ) by 


U 

s 



n 

s 


(33) 


or 


U = 
s 


[1 



+ (r /r ) 2 ] % 

s 0 s 


(34) 


The vector component of the fluid velocity normal to and measured with 
respect to the moving shock is given by 


v i * [< %. • V-vK (35) 


where q = v i + u i n . When Eqs. (32) and (34) are substituted into 

co co r co y 

Eq. (35), the following expression is obtained for the magnitude of 


V 1 = 


r - v + u (r /r ) 

S t oo oo Sq S 

[l + (r /r ) 2 ]^ 


(36) 


from which 


r 


can be obtained as 
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r = V. Ll + 
s t 1 


(r /r 

S Q S 


+ V 

co 


(r 


/r ) 
s 0 s 


(37) 


The above equations are used in the "shock- fitting" method in the 
following manner. First, the shock wave radial distance is evaluated by 
use of the Euler predictor 


n+1 

r = r" +Atr" (2 ^ k ^ NK-1) (38) 

S, S- s 

k k t, 

k 


using Eq. 
Eq. (37) 


(37) to evaluate r° . The derivative r 11 
s »- s fi 

C k 0 k 


is evaluated using the second-order central 


which appears in 


difference formula 


n 

r 

s 


(r 


k+1 


r° )/2A0 
S k-1 . 


(39) 


For the grid points immediately above (k = KS) and below (k = KS - 1) the 
shock impingement point, it is not acceptable to use Eq. (39) because 
of the discontinuity in the shock slope at the impingement point. Instead, 
the derivatives at k = KS and k = KS - 1 are evaluated using the second- 
order expressions 


n 

r 

s 





+ r 


KS + 1 



n 

r 

s 


8 


KS-1 


4r 


n 


- 3r 


n 


- r 


/3A8 


KS-1 


KS-2 


( 40 ) 
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where the radial distance to the shock impingement point (r ) is 

a 

determined in the predictor step using a second-order extrapolation 
given by 


r n = (l5 r n - 10 r n +3 r n V 8 

S a \ S KS S KS+1 S KS+2 / 


(41) 


Eqs. (40) and (41) were derived by assuming that the impingement point 
(a) is halfway between the grid points KS-1 and KS. 

The next step is to compute the pressures immediately behind 

the bow shock using the modified MacCormack scheme 


C - • f (*?.« - o - £ (<$.„ - O) - * 

(42) 


Once the pressures behind the bow shock are computed, and p? + } 

can be computed using the normal shock relations 

% 


n+1 


V, 


> ' CO 


n+1 

+ y - i 

3 Y + 1 


(43) 



(44) 


The components of the fluid velocity behind the bow shock are evaluated 
using relative velocity equations which give 
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calculated using the state relationship T = T(e, p). This completes the 
predictor step. The corrector step is identical to the predictor step 
except that the shock wave radial distance is evaluated using the modified 
Euler corrector 



and Eqs. (41) and (42) are replaced by 



( 48 ) 
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In addition, the "predicted" variables (n+1) in Eqs. (43), (44), and (45) 
are replaced by "corrected" variables (n+1) during the corrector step. 

The calculation of the boundary conditions using the "shock-fitting" 
method described above is performed before the predictor or corrector 
steps are initiated at interior grid points. All other boundary conditions 
are calculated after the predictor or corrector step is completed at 
interior grid points. The flow conditions along the outflow boundary in 
Fig. 7 are determined using a second-order extrapolation of interior grid 
point data. For example, the pressure is obtained from 


3 j,NK 3p j,NK-l ” 3p j ,NK-2 + P j,NK-3 


(49) 


Along the body surface, the following conditions are imposed 
P k,NJ = P k,NJ-l 


T = T 
k,NJ w 


^kjNJ P k,NJ /RT k,NJ 


(50.) 


U k,NJ 


= 0 


V k,NJ ° 

During some early computations, a different set of boundary conditions 
was used at the body surface. In these computations the density at the 
body surface was determined using a second-order extrapolation normal to 
the body and the pressure was found from the equation of state. These 
boundary conditions proved to be unstable, however, when applied to cases 
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where the body was highly cooled, while Eq. (50) gave stable results. 

The flow conditions along the lower boundary in Fig. 7 are determined 
using either simple reflection about an axis of symmetry which is placed 
midway between the k = 1 and k = 2 rows of grid points or a second-order 
extrapolation of interior grid point data for the case of supersonic 
outflow conditions. 


Initial Conditions 

In all computations performed thus far, the blunt body flow without 
the impinging shock was computed first. The initial conditions for this 
calculation are obtained by using an approximate curve fit for the loca- 
tion and shape of the bow shock along with a Newtonian pressure distribu- 

12 

tion at the body. The approximate curve fit of Billig is used to find 
r and r along the shock. Eqs. (36), (43), (44), and (45) are then used 

S Sq 

(with r set equal to zero) to find the initial conditions immediately 
S t 

behind the shock. The initial flow conditions at the wall are obtained 
using the known wall temperature in conjunction with the pressures from 
the Newtonian pressure distribution. The densities at the wall are obtained 
from the equation of state and the velocities are set equal to zero. The 
initial flow conditions at interior grid points are obtained by assuming 
a linear variation between the flow conditions immediately behind the 
bow shock and the wall conditions. 

After a "steady-state" solution is achieved for the undisturbed 
blunt body flow, the impinging shock is introduced by resetting the free- 
stream flow variables above the impingement point equal to the values 
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which exist behind the desired oblique impinging shock. The computation 
is then carried to a new "steady-state". 

Stability 

Since the present computational method is explicit, the maximum time 
increment (At) must be limited to ensure stability. , For high Reynolds 
number flows, the maximum time increment is determined from the usual 

g 

C.F.L. condition . For low Reynolds number flows, the following equation 
13 

derived by Li can be used: 

At S P(AT P 2 


8p, 


(51) 
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NUMERICAL RESULTS 


The computational methods developed in this study have been applied 
to several 2-D viscous, blunt body flows with an impinging shock. Both 
the "shock-capturing" and "shock-fitting" methods have been used. 


Case I. ("Shock-Capturing" method) 


The first test case consists of nitrogen flowing over a circular 
cylinder at low Reynolds number. The flow conditions for this case 
are 


M = 4.2 

CD 

Re D = 200 

CO 


2 

p = .0735 newtons/m 

oo 

T = 96.1 °K 

CO 


(52) 


Pr = .687 


D = .6096 m 


y = 1.4 


T = T = 435.2 °K 
w t 

OO 


The "shock-capturing" method is ideally suited for this computation since 
the bow shock is relatively thick for low values of Reynolds number. For 
this case, the outer computational boundary was a circle with a radius 
equal to the diameter of the cylinder. The lower boundary (k=l) was 
located along a ray 1° below the undisturbed stagnation streamline and 
the outflow boundary (k=NK) was located along a ray 39.1° above the 
undisturbed stagnation streamline. A mesh consisting of 21 equally 
spaced grid points in both the y and z directions was used in this 


computation . 
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Initially, the blunt body flow without the impinging shock was 
computed using freestream initial conditions at all grid points. This 
computation took 1425 time steps to reach a "steady~state" and used 
approximately 24.2 minutes of CPU time on the IBM 360-65 computer. The 
results of this computation are seen in the contour plots of Figs. 9, 10, 
and 11 which show lines of constant Mesh number, temperature, and pressure, 
respectively. The finite thickness of the bow shock is evident in these 
figures. 

The variation of the computed temperature along the stagnation 
streamline is shown in Fig. 12, along with the experimentally measured 
wire temperatures of Ref. 14. Although the wire temperatures cannot 
be compared directly with the computed static temperatures, they can be 
used to indicate the location of the bow shock. The wire temperatures 
have been normalized so that their values in the freestream and at the 
body agree with the computed static pressures at these points. The 
location of the bow shock in the present computation falls between the 
experimentally determined locations in Fig. 12. This is the desired 
result, since the present flow conditions are between the flow condi- 
tions of the two experimental tests of Ref. 14. 

An impinging shock was then introduced into the "steady-state" blunt 
body computation along a ray 20° above the stagnation streamline of the 
undisturbed flow field. This was accomplished by resetting the freestream 
flow variables, along the outer computational boundary above the impinge- 
ment point, equal to the values which exist behind a 20° shock. The 
ratio of the pressures across this impinging shock was p^ /p^ = 2.24. 

CD CO 




Fig. 9. Constant Mac 
Case I (no 
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Constant pressure lines, 
Case I (no impingement). 


TEMPERATURE, T A 
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TRANSFORMED COORDINATE, z 


Fig. 12. Variation of temperature along 
stagnation streamline. 
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The computation was then carried to a new "steady-state" which took an 
additional 2325 steps and used approximately 39.5 minutes of CPU time 
on the IBM 360-65 computer. The results of this computation are shown 
in the contour plots of Figs. 13, 14, and 15. The large effect of the 



Fig. 13. Constant Mach number lines. Case I 
(20° shock impingement). 


Fig. 14. Constant temperature lines. Case I . 

(20° shock impingement) . 

incident shock on the undisturbed flow field is evident when these 
figures are compared with the previous contour plots. Comparisons of 
the surface pressures and heat transfers between the undisturbed 


computation and the present shock impingement computation are shown 
in Figs. 16 and 17. The incident shock causes an increase in wall 
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Fig. 15. Constant pressure lines. Case I . 

(20° shock impingement). 

pressure of 72% and an increase in heat transfer of 49%, over the un- 
disturbed values at the point where the outflow boundary (k=21) inter- 
sects the body. 


WALL PRESSURE 
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Fig. 16. Comparison of wall pressures. 



HEAT TRANSFER 
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Fig. 17. Comparison of heat transfers. 
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Case II. ("Shock-Fitting" Method) 

This test case consists of air flowing over a circular cylinder with 


the following conditions 


M = 12.145 

CD 

p =19.8 newtons/m 

r CD 

Re D = 19,838 

CO 

T = 254 °K 

CO 

Pr = .72 

D = .3048 m 

\ = 1.4 

T = 1445 °K 
w 


These flow conditions are identical to those that will exist normal to 
the leading edge of the Orbiter wing (sweep angle of 45°) at 61 km 
(200,000 ft). 

The "shock-fitting" method was used for this computation since the 
moderately high Reynolds number causes the bow shock to be relatively 
thin. A mesh consisting of 51 equally spaced grid points in the y 
direction and 31 equally spaced grid points in the z direction was used. 

The lower outflow boundary (k=l) was located along a ray 61.5° below the 
undisturbed stagnation streamline while the upper outflow boundary was 
located along a ray 88.5° above the undisturbed stagnation streamline. 

These boundaries were placed so that the flow passing out of each of 
them was supersonic everywhere except in the boundary layer. 

As before, the blunt body flow without the impinging shock was 
computed initially. The results of this computation for 8=0 are shown in the 
contour plots of Figs. 18, 19, and 20. Only the results from k=21 to k=51 
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are plotted. The stagnation streamline is located halfway between k=21 
and k=22. The boundary layer is clearly evident in these figures. This 

undisturbed computation was repeated using both the exponential stretching 

/ 

(Eq. (16), with 3=3) and the logarithmic stretching (Eq. (Al), with 
3' = 1.2 and a = 0) . A comparison of the temperature profiles at the 
stagnation point for the three different computations is shown in Fig. 21. 
Included in this figure is the temperature profile computed by the 
BLIMP^ boundary layer program, assuming a perfect gas. Excellent agree- 
ment is achieved for all computations. Comparisons of the temperature 
and density profiles at various points around the cylinder for the 
logarithmic stretching case are shown in Figs. 22 and 23. Again, the 
agreement between the present results and the BLIMP program results is 
very good. 

An impinging shock wave was then introduced into the previous "steady 

state" undisturbed computation. This relatively weak shock wave was 

inclined to the freestream at an angle of 5%° and impinged on the bow 

shock at 0 = 39° (halfway between k=34 and k=35 along j=l). The ratio 

of the pressures across this impinging shock was / p^ = 1.42. The 

00 00 

computation was then carried to a new "steady-state" and the resulting 
contour plots for 3=0 are shown in Figs. 24, 25, and 26. These 
figures clearly show the shear layer which originates from the impingement 
point and passes out the upper outflow boundary. This shock interaction 
pattern corresponds to Edney's Type III. The effect of the impinging 
shock on the shape of the bow shock is shown in Fig. 27. 

When the present shock impingement case was recomputed using 
exponential stretching ( 3 = 3), an instability developed which caused 
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Fig. 21. Temperature profiles at stagnation point. 





LOGARITHMIC STRETCHING (p* 
BLIMP PROGRAM 


o 

it 
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Fig. 23. Density profiles 
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Fig. 27. Shock shapes 
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the solution to "blow-up". This instability was believed to be caused 
by the large grid spacing near the bow shock for the exponential type of 
stretching. This difficulty was circumvented by using logarithmic 
stretching with P' = 1.12 and a. = 0. This value of (3' gives the same 
refinement at the body as the exponential type of stretching with (3 = 3, 
but the grid spacing near the bow shock is smaller for the logarithmic 
type of stretching. 

Comparisons of the surface pressures and heat transfers between the 
undisturbed computation and the present shock impingement computation 
with logarithmic stretching are shown in Figs. 28 and 29. The incident 
shock causes an increase in wall pressure of 33% over the undisturbed 
values at the point where the outflow boundary (k = 51) instersects the 
body. The large singular value of heat transfer near 0=5° for the 
shock impingement case is the result of a numerical wave which emanates 
from the impingement point and hits the body. This numerical wave is 
evident in Fig. 25 and is believed to be caused by the finite values of 
shock velocity which remain near the impingement point even after the 
"steady-state" solution is reached. These finite values of shock 
velocity alternate signs between the predictor and corrector steps but 
have the same magnitude so that the shock position does not change after 
a complete time step is computed. Work is underway to eliminate this 
so-called "chattering" effect. 

Attempts to introduce a stronger impinging shock or attempts to move 
the impingement point further into the subsonic region for the present flow 
conditions have been unsuccessful to date. Work is continuing in an 


effort to overcome these difficulties. 



49 


o NO SHOCK IMPINGEMENT 
A 5 1/2 ° SHOCK IMPINGEMENT AT e = 39° 



-60 -40 -20 0 20 40 60 80 

ANGLE, 9 , DEG. 


Fig. 28. Comparison of wall pressures. 
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ONO SMOCK IMPINGEMENT 


A 5 1/5° SHOCK IMPINGEMENT AT 0= 39° 
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Fig. 29. Comparison of heat transfers. 
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APPENDIX A 


The transformation given below is similar to one developed by Roberts 
and has been found in the present study to give better computational results 
than the previous transformation of Li^ when an impinging shock is present. 
This is due to the fact that the grid spacing near the bow shock using 
this transformation is smaller than the grid spacing of Li's transformation 
for an equal grid spacing at the body. The new transformation can be 
used to refine the mesh near the body (Oi = 0), as was done in this study, 
or it can be used to refine the mesh at both the body and bow shock 
(a = %) . The equations for this transformation are 



z 


a + (l - a) 



B' + z (2a + 1) - 2 cA 

B' - z (2a + l) + 2a J 



(Al) 


t = t 


The relations between the partial derivatives are given by 


d_ 

dy 


d_ 

% 


(A2) 


d_ 

dz 


[b' + z(2a + 1) - 20t] [b’ - z(2q + 1) + 2 a] ln^ ~gT~7~Y^ 


2 p' (i - a) (2a + l) 



d_ 

dz 



54 


where z - g 


O' +2 0) ( 

+ i.y ■ “ g> 

f 2a 

-i) e 

(2a + l) 

z - a 


.-(f-H) 1 '". 


(A3) 


and j3' is related to the thickness (c) of the "boundary layer" where 
clustering is desired by 


0' - (1 - c) 


0 < c < 1 


(A4) 



